##########################################
########################################## 
## STUDY 3
##########################################
########################################## 
setwd(getwd())
source("libraries.R")

# load data for study
data_raw <- haven::read_dta(here("Data", "study3.dta"))

##########################################
########################################## 
## Recode Variables
##########################################
########################################## 
df <- data_raw %>% 
  mutate(female = case_when(
    gender == 1 ~ 0,
    gender == 2 ~ 1,
    TRUE ~ NA_real_
  ), 
  race = case_when(
    race == 1 ~ 1,
    TRUE ~ 0
  )) %>% 
  rename(edu = educ,
         party_ID = Q6a) %>% 
  
  #Party Identification (1 = Democrat; 0 = Republican)
  mutate(repdem = case_when(
    party_ID == 1 | party_ID == 2 | party_ID == 3 ~ 1,
    party_ID == 5 | party_ID == 6 | party_ID == 7 ~ 0,
    TRUE ~ NA_real_
  ),
  
  #Dummy indicating whether participant is Independent or not
  independents = case_when(
    party_ID == 4 ~ 1,
    TRUE ~ NA_real_
  )) %>% 
  
  #Code "Don't Know" as missing
  mutate_at(vars(Q7_1:Q7_4, Q8_1:Q8_8, Q10_1:Q11_7, Q13_1:Q20_2), ~na_if(., 977))  %>% 
  
  #Measure of straight-lining
  mutate(straight_line = case_when(
    Q7_1 == 1 & Q7_2 == 1 & Q7_3 == 1 & Q7_4 == 1 ~ 1,
    Q7_1 == 2 & Q7_2 == 2 & Q7_3 == 2 & Q7_4 == 2 ~ 1,
    Q7_1 == 3 & Q7_2 == 3 & Q7_3 == 3 & Q7_4 == 3 ~ 1,
    Q7_1 == 4 & Q7_2 == 4 & Q7_3 == 4 & Q7_4 == 4 ~ 1,
    Q7_1 == 5 & Q7_2 == 5 & Q7_3 == 5 & Q7_4 == 5 ~ 1,
    TRUE ~ 0
  )) %>% 

  #Need for Chaos
  mutate(need_chaos = rowMeans(select(., Q8_1:Q8_7)), 
                        na.rm = TRUE) %>% 
  
  #Political Violence Intentions
  mutate(violentact = rowMeans(select(., c("Q10_6", "Q10_7", "Q10_8", "Q11_1", "Q11_2", "Q11_3", "Q11_4",
                                           "Q11_5", "Q11_6", "Q11_7"))), 
                               na.rm = TRUE) %>% 
  
  #Political Rumors: Believing and Sharing of DEMocratic and REPublican rumors
  mutate(true_dem = rowMeans(select(., c("Q13_1","Q14_1","Q15_1"))), na.rm = TRUE) %>% 
  mutate(share_dem = rowMeans(select(., c("Q13_2","Q14_2","Q15_2"))), na.rm = TRUE) %>% 
  mutate(true_rep = rowMeans(select(., c("Q18_1","Q19_1","Q20_1"))), na.rm = TRUE) %>% 
  mutate(share_rep = rowMeans(select(., c("Q18_2","Q19_2","Q20_2"))), na.rm = TRUE) %>%  
  
  #Standardize variables
  mutate(across(.cols = c(need_chaos, violentact, age, edu, true_dem, share_dem, true_rep, share_rep),
                ~ (scale(.) %>%  as.vector),
                .names = "{.col}_zscore")) %>% 
  
  #Scale variables to range from 0 to 1
  mutate(across(.cols = c(need_chaos, violentact, age, edu, true_dem, share_dem, true_rep, share_rep),
                ~ (zero1(.) %>% as.vector),
                .names = "{.col}01"))

##########################################
########################################## 
## Sample Description 
## (See Supplementary Materials SM1, Table S2)
##########################################
########################################## 
#Gender
print(paste(round(mean(df$female, na.rm = T), 2)*100 , "percent of participants are female"))
#Age
print(paste("Mean age of participants:" , round(mean(df$age), 2), 
            ", SD:", round(sd(df$age) , 2) ))
#Education
print(paste("Less than high school:", round(prop.table(table(df$edu)), 2)[1]*100, "percent"))
print(paste("High school graduate:", round(prop.table(table(df$edu)), 2)[2]*100, "percent"))
print(paste(" Some college, but no degree:", round(prop.table(table(df$edu)), 2)[3]*100, "percent"))
print(paste("2 year college degree:", round(prop.table(table(df$edu)), 2)[4]*100, "percent"))
print(paste("4 year college degree:", round(prop.table(table(df$edu)), 2)[5]*100, "percent"))
print(paste("Post-grad:", round(prop.table(table(df$edu)), 2)[6]*100, "percent"))
#Ethnicity
print(paste("White/caucasian:", round(prop.table(table(df$race)), 2)[2]*100, "percent"))

##########################################
########################################## 
## Factor Analysis, etc.: Need for Chaos 
## (See Supplementary Materials SM2, Table S3)
##########################################
########################################## 
#Cronbach's alpha
nfc <- df %>% 
  select(Q8_1:Q8_7) %>% 
  psych::alpha(.) 
print(paste("Cronbach's alpha, Need for Chaos:" , round(nfc$total$raw_alpha,2)))  

#Confirmatory factor analysis
cfa <- df %>% 
  select(Q8_1:Q8_7) %>% 
  na.omit()

#One-factor model (See Table S7)
m1  <- ' nfc_oneDim  =~ Q8_1 + Q8_2 + Q8_3 + Q8_4 + Q8_5 + Q8_6 + Q8_7  
        Q8_3~~Q8_4
        Q8_3~~Q8_5
        Q8_4~~Q8_5'
m1 <- lavaan::cfa(m1, data=cfa) 
m1_fit <- summary(m1 , fit.measures=TRUE) 

print(paste("One-factor model, Chi-square:" , round(m1_fit$FIT["chisq"],2 )))
print(paste("One-factor model, CFI:" , round(m1_fit$FIT["cfi"],2 )))
print(paste("One-factor model, TLI:" , round(m1_fit$FIT["tli"],2 )))
print(paste("One-factor model, RMSEA:" , round(m1_fit$FIT["rmsea"],2 )))

#Two-factor model
m2  <- ' nfc_apol  =~ Q8_1 + Q8_2 + Q8_6 + Q8_7  
         nfc_pol  =~ Q8_3 + Q8_4 + Q8_5'
m2 <- lavaan::cfa(m2, data=cfa) 
m2_fit <- summary(m2 , fit.measures=TRUE) 

print(paste("Two-factor model, Chi-square:" , round(m2_fit$FIT["chisq"],2 )))
print(paste("Two-factor model, CFI:" , round(m2_fit$FIT["cfi"],2 )))
print(paste("Two-factor model, TLI:" , round(m2_fit$FIT["tli"],2 )))
print(paste("Two-factor model, RMSEA:" , round(m2_fit$FIT["rmsea"],2 )))

#Compare one- vs two-factor model
compare <- anova(m1, m2, method = "satorra.bentler.2010")

print(paste("Chi-square difference:" , round(compare[2,5],2)))
print(paste("Chi-square difference, p-value:" , compare[2,7]))

##########################################
########################################## 
## Remove 5% with highest Need for Chaos? (SM8)
## When set to TRUE, used for producing Fig S8, Fig S9
##########################################
##########################################
# Set to TRUE for excluding top 5% of participants with highest Need for Chaos
exclude_obs <- FALSE

if (exclude_obs == TRUE) {
  
  df <- df %>% 
    filter(need_chaos_zscore < quantile(need_chaos_zscore, probs = .95, na.rm = TRUE))
  
  print("Excluding top 5% observations with highest levels of the Need for Chaos")
  
}

##########################################
########################################## 
## Analysis - Fig. 3 of Main Text
## Note (Q24 == 1 excludes trolls while straight_line == 0 excludes straight-liners)
##########################################
########################################## 
m_share_dem <- lm(share_dem01 ~ need_chaos_zscore + repdem + female + edu_zscore + age_zscore + race , data = df[df$Q24 == 1 & df$straight_line == 0, ] )
share_dem <- tidy(m_share_dem) %>% 
  filter(term == "need_chaos_zscore" | term == "repdem") %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_dem",
         study = "study3")

m_share_rep <- lm(share_rep01 ~ need_chaos_zscore + repdem +  female + edu_zscore + age_zscore + race , data = df[df$Q24 == 1 & df$straight_line == 0, ] )
share_rep <- tidy(m_share_rep) %>% 
  filter(term == "need_chaos_zscore" | term == "repdem") %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_rep",
         study = "study3")

m_true_dem <- lm(true_dem01 ~ need_chaos_zscore + repdem  +  female + edu_zscore + age_zscore + race , data = df[df$Q24 == 1 & df$straight_line == 0, ] )
true_dem <- tidy(m_true_dem) %>% 
  filter(term == "need_chaos_zscore" | term == "repdem") %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_dem",
         study = "study3")

m_true_rep <- lm(true_rep01 ~ need_chaos_zscore + repdem  +  female + edu_zscore + age_zscore + race , data = df[df$Q24 == 1 & df$straight_line == 0, ] )
true_rep <- tidy(m_true_rep) %>% 
  filter(term == "need_chaos_zscore" | term == "repdem") %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_rep",
         study = "study3")

study3_coef <- rbind(share_dem, share_rep, true_dem, true_rep)

#save coefficients for figure 
save(study3_coef,file=here("Data_recoded","study3_fig3.R"))

#regression table for Supplemental Materials (Table S10)
stargazer(m_share_dem , m_share_rep , m_true_dem , m_true_rep , 
          align=TRUE,
          no.space=TRUE,
          column.labels=c("Share - Dem", "Share - Rep" , "True - Dem" , "Share - Rep"),
          covariate.labels = c("Need for Chaos (std.)" , "Party ID (1 = Dem.)" , "Women" , "Education" , "Age  (std.)" , "Ethnicity (1 = non-white)" ,  "Intercept"),
          keep.stat = c("n","ser"),
          dep.var.labels.include = FALSE)

##########################################
########################################## 
## Analysis - Fig. 4 of Main Text
##########################################
########################################## 
# Fig 4: Democratic Voters  
df_dem <- df %>% 
  filter(repdem == 1 & Q24 == 1 & straight_line == 0 )

m_share_dem <- lm(share_dem01 ~ need_chaos_zscore +  female + edu_zscore + age_zscore + race , data = df_dem )
share_dem <- tidy(m_share_dem) %>% 
  filter(term == "need_chaos_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_dem",
         study = "study3")

m_share_rep <- lm(share_rep01 ~ need_chaos_zscore +  female + edu_zscore + age_zscore + race , data = df_dem )
share_rep <- tidy(m_share_rep) %>% 
  filter(term == "need_chaos_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_rep",
         study = "study3")

m_true_dem <- lm(true_dem01 ~ need_chaos_zscore +  female + edu_zscore + age_zscore + race , data = df_dem )
true_dem <- tidy(m_true_dem) %>% 
  filter(term == "need_chaos_zscore") %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_dem",
         study = "study3")

m_true_rep <- lm(true_rep01 ~ need_chaos_zscore +  female + edu_zscore + age_zscore + race , data = df_dem )
true_rep <- tidy(m_true_rep) %>% 
  filter(term == "need_chaos_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_rep",
         study = "study3")

study3_coef_dem <- rbind(share_dem, share_rep, true_dem, true_rep)

# Fig 4: Republican Voters 
df_rep <- df %>% 
  filter(repdem == 0 & Q24 == 1 & straight_line == 0 )

m_share_dem <- lm(share_dem01 ~ need_chaos_zscore +  female + edu_zscore + age_zscore + race , data = df_rep )
share_dem <- tidy(m_share_dem) %>% 
  filter(term == "need_chaos_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_dem",
         study = "study3")

m_share_rep <- lm(share_rep01 ~ need_chaos_zscore +  female + edu_zscore + age_zscore + race , data = df_rep )
share_rep <- tidy(m_share_rep) %>% 
  filter(term == "need_chaos_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_rep",
         study = "study3")

m_true_dem <- lm(true_dem01 ~ need_chaos_zscore +  female + edu_zscore + age_zscore + race , data = df_rep )
true_dem <- tidy(m_true_dem) %>% 
  filter(term == "need_chaos_zscore") %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_dem",
         study = "study3")

m_true_rep <- lm(true_rep01 ~ need_chaos_zscore +  female + edu_zscore + age_zscore + race , data = df_rep )
true_rep <- tidy(m_true_rep) %>% 
  filter(term == "need_chaos_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_rep",
         study = "study3")

study3_coef_rep <- rbind(share_dem, share_rep, true_dem, true_rep)

#save coefficients for figure 
save(study3_coef_dem,file=here("Data_recoded","study3_fig4_dem.R"))
save(study3_coef_rep,file=here("Data_recoded","study3_fig4_rep.R"))

##########################################
########################################## 
## Analysis - Fig. S2 of Supplementary Materials
##########################################
########################################## 
df_independents <- df %>% 
  filter(independents == 1 & Q24 == 1 & straight_line == 0 )

m_share_dem <- lm(share_dem01 ~ need_chaos_zscore +  female + edu_zscore + age_zscore + race , data = df_independents )
share_dem <- tidy(m_share_dem) %>% 
  filter(term == "need_chaos_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_dem",
         study = "study3")

m_share_rep <- lm(share_rep01 ~ need_chaos_zscore +  female + edu_zscore + age_zscore + race , data = df_independents )
share_rep <- tidy(m_share_rep) %>% 
  filter(term == "need_chaos_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "share_rep",
         study = "study3")

m_true_dem <- lm(true_dem01 ~ need_chaos_zscore +  female + edu_zscore + age_zscore + race , data = df_independents )
true_dem <- tidy(m_true_dem) %>% 
  filter(term == "need_chaos_zscore") %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_dem",
         study = "study3")

m_true_rep <- lm(true_rep01 ~ need_chaos_zscore +  female + edu_zscore + age_zscore + race , data = df_independents )
true_rep <- tidy(m_true_rep) %>% 
  filter(term == "need_chaos_zscore" ) %>% 
  select("term","estimate", "std.error") %>% 
  mutate(dv = "true_rep",
         study = "study3")

study3_coef_independents <- rbind(share_dem, share_rep, true_dem, true_rep)

#save coefficients for figure 
save(study3_coef_independents,file=here("Data_recoded","study3_figS2.R"))

##########################################
########################################## 
## Analysis - Supplementary Materials SM8
##########################################
########################################## 
#Does Need for Chaos predict intentions to engage in political violence? 
m1 <- lm(violentact01 ~ need_chaos_zscore + repdem  +  female + edu_zscore + age_zscore + race , data = df[df$Q24 == 1 & df$straight_line == 0, ])

#regression table for supplemental materials (Table S28)
stargazer(m1 , 
          align=TRUE,
          no.space=TRUE,
          column.labels=c("Political Violence Intentions"),
          covariate.labels = c("Need for Chaos (std.)" , "Party ID (1 = Dem.)" , "Women" , "Education" , "Age  (std.)" , "Ethnicity (1 = non-white)" ,  "Intercept"),
          keep.stat = c("n","ser"),
          dep.var.labels.include = FALSE)

########################
########################
########################
## Panel Study of Wave 1,2,3 (see Supplementary Materials SM3)
########################
########################
########################

##########################################
########################################## 
## Recode Variables
##########################################
########################################## 
df_panel <- df %>% 
  mutate_at(vars(Q1_1_2:Q1_8_2, Q1_1_3:Q1_8_3, Q3_1_2:Q8_2_2), ~na_if(., 977)) %>% 
  
  #Need for Chaos in Wave 1, 2, 3
  mutate(nfc1 = rowMeans(select(., c("Q8_1", "Q8_2", "Q8_3", 
                                           "Q8_4", "Q8_5", "Q8_6", "Q8_7"))), 
                               na.rm = TRUE) %>% 
  mutate(nfc2 = rowMeans(select(., c("Q1_1_2", "Q1_2_2", "Q1_3_2", "Q1_4_2", 
                                     "Q1_5_2", "Q1_6_2", "Q1_7_2"))), 
                         na.rm = TRUE) %>% 
  mutate(nfc3 = rowMeans(select(., c("Q1_1_3", "Q1_2_3", "Q1_3_3", 
                                     "Q1_4_3", "Q1_5_3", "Q1_6_3", "Q1_7_3"))), 
                         na.rm = TRUE) %>% 
  
  ##Political Rumors: Believing and Sharing in Wave 2
  mutate(true_2dem = rowMeans(select(., c("Q3_1_2","Q4_1_2","Q5_1_2"))), na.rm = TRUE) %>% 
  mutate(share_2dem = rowMeans(select(., c("Q3_2_2","Q4_2_2","Q5_2_2"))), na.rm = TRUE) %>% 
  mutate(true_2rep = rowMeans(select(., c("Q6_1_2","Q7_1_2","Q8_1_2"))), na.rm = TRUE) %>% 
  mutate(share_2rep = rowMeans(select(., c("Q6_2_2","Q7_2_2","Q8_2_2"))), na.rm = TRUE) %>% 
  
  ##Changes in NFC and Rumor Endorsement between Wave 1 and Wave 2, standardization and 0-1 coding
  mutate(nfc_change = nfc2 - nfc1,
         share_dem_change = share_2dem - share_dem,
         share_rep_change = share_2rep - share_rep,
         true_dem_change = true_2dem - true_dem,         
         true_rep_change = true_2rep - true_rep,         
         
         nfc_change_std = standardize(nfc_change),
         share_dem_change_01 = zero1(share_dem_change),
         share_rep_change_01 = zero1(share_rep_change),
         true_dem_change_01 = zero1(true_dem_change),
         true_rep_change_01 = zero1(true_rep_change),
         )

##########################################
########################################## 
#Table S14: Did changes in Need for Chaos 
##predict changes in rumor endorsement from 
##survey wave 1 to survey wave 2?
##########################################
########################################## 
reg_dem <- lm(share_dem_change_01 ~ nfc_change_std, data = df_panel)
reg_rep <- lm(share_rep_change_01 ~ nfc_change_std, data = df_panel)
reg_dem2 <- lm(true_dem_change_01 ~ nfc_change_std, data = df_panel)
reg_rep2 <- lm(true_rep_change_01 ~ nfc_change_std, data = df_panel)

stargazer::stargazer(reg_dem,reg_rep,reg_dem2,reg_rep2,
                     align=TRUE,
                     no.space=TRUE,
                     keep.stat = c("n"),
                     dep.var.labels.include = TRUE,
                     covariate.labels = c("Change in NFC (std.)", "Intercept"))

##########################################
########################################## 
#Analysis: Stability of NFC across survey waves
##########################################
########################################## 
#Correlation Measures
df_panel <- df_panel %>% 
  mutate(across(.cols = c(nfc1, nfc2, nfc3),
                ~ (zero1(.) %>% as.vector),
                .names = "{.col}_01")) %>% 
    select(RecordNo, nfc1_01 , nfc2_01 , nfc3_01)

df_long <- pivot_longer(data = df_panel , cols = nfc1_01:nfc3_01 , names_to = "wave")

df_panel <- df_panel[complete.cases(df_panel),]
print(paste("Need for Chaos correlation between wave 1 & 2:" , round(cor(df_panel$nfc1_01 , df_panel$nfc2_01),2)))
print(paste("Need for Chaos correlation between wave 1 & 3:" , round(cor(df_panel$nfc1_01 , df_panel$nfc3_01),2)))

#Figure S3: Need for Chaos across three survey waves
w1 <- df_long %>% filter(.,wave == "nfc1_01")
q_w1 <- quantile(w1$value , c(.99 , .90, .25) , na.rm = TRUE)

w2 <- df_long %>% filter(.,wave == "nfc2_01")
q_w2 <- quantile(w2$value , c(.99 , .90, .25) , na.rm = TRUE)

w3 <- df_long %>% filter(.,wave == "nfc3_01")
q_w3 <- quantile(w3$value , c(.99 , .90, .25) , na.rm = TRUE)

top_1 <- data.frame(mean = c(q_w1[1]  , q_w2[1]  , q_w3[1]), wave = c("nfc1_01","nfc2_01","nfc3_01"))
top_10 <- data.frame(mean = c(q_w1[2]  , q_w2[2]   , q_w3[2]), wave = c("nfc1_01","nfc2_01","nfc3_01"))
bottom_25 <- data.frame(mean = c(q_w1[3]  , q_w2[3]   , q_w3[3]), wave = c("nfc1_01","nfc2_01","nfc3_01"))

figure <- ggplot() + 
  geom_line(data=df_long, aes(x=as.factor(wave), y=value , group = RecordNo), colour="black" , size = 1.4 , alpha=0.02 , position=position_jitter(w=0.02, h=0.03)) + 
  geom_line(data=top_1, aes(x=as.factor(wave), y=mean , group = 1), color='black' , size = 1.6) +
  geom_line(data=top_10, aes(x=as.factor(wave), y=mean , group = 1), color='black' , size = 1.6) +
  geom_line(data=bottom_25, aes(x=as.factor(wave), y=mean , group = 1), color='black' , size = 1.6) +
  annotate("text", x = .77, y = q_w1[1], label = "W1: top1%" , size = 3.2) +
  annotate("text", x = .77, y = q_w1[2], label = "W1: top10%" , size = 3.2) +
  annotate("text", x = .75, y = q_w1[3], label = "W1: bottom25%" , size = 3.2) +
  xlab("") + ylab("Need for Chaos") + labs(title = " "  ) +
  scale_x_discrete(breaks=c("nfc1_01","nfc2_01","nfc3_01"),
                   labels=c("Wave 1", "Wave 2 (+1-2 months)", "Wave 3 (+7-10 months)")) 

ggsave(here::here("Figures", "figS3.pdf"),
       width = 8, height = 6)
